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different transients apart. Two algorithms are evaluated. One method is based on higher 
order statistics while the other is based on signal subspace techniques. Extensive 
simulations using synthetic transients were conducted to establish the performance of 
each algorithm in terms of discrimination and TDOA estimation. It was found that the 
bispectral algorithm gave better TDOA estimation at low SNRs while the subspace 
algorithm gave better TDOA estimation at high SNRs. For discrimination, it was found 
that the subspace algorithm gave consistant false alarm rates at all SNRs while the false 
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EXECUTIVE SUMMARY 


In this thesis we have developed and compared two algorithms, namely the 
bispectrum and subspace linear phase detectors. These algorithms were developed for 
the purposes of transient discrimination and Time-Difference-Of-Arrival (TDOA) 
estimation. They are to be used as part of a transient tool suite to aid in the estimation of 
a submarine’s position. Two performance measures were used to evaluate the 
algorithms, namely the probability of correct TDOA (Pr) and the probability of correct 
discrimination (Pp;). 

In general, it can be said that for TDOA, the bispectral linear phase detector gave 
better results at low SNRs while the subspace linear phase detector worked better at the 
higher SNRs. 

For discrimination, 1t was found that the bispectral discriminator gave higher Pp; 
than the subspace discriminator. However, the probability of false discrimination (Pg) of 
the bispectral discriminator increased at higher SNRs while the subspace discriminator 
gave a constant Pr, for all SNRs evaluated. For discrimination, the advantage of having a 
constant Pg is desirable. Therefore the subspace discriminator is the best option even 
although it produced lower Pp; than the bispectral discriminator. It was also found that 
there are design trade-offs between processing speed and performance that need to be 
made. For the bispectral linear phase detector, this trade-off is in terms of threshold gain; 


for the subspace linear phase detector, this trade-off is in terms of correlation matrix size. 
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I. INTRODUCTION 


The work of this thesis forms part of the ongoing effort to automate the detection, 
discrimination and. Time-Difference-Of-Arrival (TDOA) estimation of transient signals. 
At present, transient signals are located and processed manually, making it a very time 
consuming procedure. It is therefore more desirable to have an automated system that 
can deliver the same or better results than the present manual approach. 

For this thesis, it 1s assumed that the detection of the transient has already taken 
place; therefore, the objective of the thesis 1s to look at discrimination and TDOA 
estimation only. Two algorithms, namely the bispectral linear phase detector and the 


subspace linear phase detector, are developed for this purpose. 
A. THE OPERATIONAL MISSION 


A typical mission considered for this thesis starts off with an anti-submarine 
warfare (ASW) aircraft laying a sonabuoy field in the vicinity or on the predicted path of 
a submarine. Usually each sonabuoy in the field consists of calibrated omni-directional 
passive acoustic sensors. The ASW aircraft collects acoustic data on a target submarine 
as It passes through the sonabuoy field. 

Figure 1 shows a typical V-shaped sonabuoy field that can be used in this type of 
mission. When this shape is used, the target must pass through the apex of the V of the 
field (as close as possible) in order to obtain the most accurate results. The data received 
at the sonabuoys is transmitted to the aircraft, which in turn records the data on tape for 


mission post processing. 





Figure 1. Typical Sonabuoy Field. 


The data collected on the tapes 1s processed at an onshore location in order to estimate the 
acoustic Sound Pressure Levels (SPL) of tonals, broadband signals and transients that are 
emitted from the submarine. 

To obtain precise SPL of a submarine, an accurate track is required so that the 
pressure levels received at the sonabuoys can be projected back to estimate the pressure 
levels at the target. For best results, the position of each sonabuoy must be known, the 
target must be on a course that passes through the middle of the field, and an accurate 


estimation of TDOA of acoustic signals must be made. 
B. DATA PROCESSING 


The data collected from the sonabuoys and stored on the tape is processed by first 
generating a baseband sonogram. А sonogram is a frequency versus time plot of the 
acoustic pressure levels at a sonabuoy. From these sonograms an analyst can identify and 
highlight contact events of interest, which may include signals, such as narrowband 
tonals and transients. The events identified are then processed to obtain an estimate of 
the target’s track. Once the best estimate of the target track has been made, the analyst 
can conduct SPL analysis on all types of previously identified contact signals. The final 
desired results are the SPL signature characteristics for a given target with respect to 


aspect angle. 


At present, transient analysis is a manually intensive operation. Typically, the 
analyst marks transient events of interest from the sonogram display of a sonabuoy and 
then searches the sonograms of the other sonabuoys for the same transient. This can be a 
very time consuming process since in a single mission there can be hundreds of transients 
within a single sonogram. These must be matched up to those on the sonograms from the 


other sonabuoys in the sonabuoy field. 
С. THESIS GOAL 


The goal of this thesis 1s to develop and evaluate two algorithms that can be used 
for transient discrimination and. Time-Difference-Of-Arrival (TDOA) estimation. The 
two algorithms are the bispectral linear phase detector and the subspace linear phase 


detector. 
D. THESIS OUTLINE 


The remainder of this thesis is organized as follows. Chapter Il discusses 
transients and transient processing and introduces the problems of interest, namely 
TDOA estimation and transient discrimination. Chapter III presents the relevant signal 
models and the analysis of the measured signals in terms of second- and third-order 
moments. This chapter details the formulation of the bispectral linear phase detector and 
discriminator. Chapter IV discusses signal subspace techniques and their application to 
the problem of TDOA estimation and transient discrimination. The outcome of this 
chapter is the formulation of the subspace TDOA estimator and discriminator. Chapter V 
examines the results achieved by both the bispectral and subspace TDOA estimators and 
discriminators. Extensive simulations are conducted using synthetic transients to produce 
performance curves and to study the utility of each technique in discrimination and 
TDOA of transients. Chapter VI provides conclusions and recommendations for future 
work based on the results presented. Appendix A and Appendix B show the 


discrimination results for the different combinations of signals not shown in Chapter V. 
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II. TRANSIENT PROCESSING 


A. TRANSIENT SIGNALS 


Transients of interest here are short duration wideband acoustic signals. 
Accordingly transients can be of varying shapes and lengths, lasting anywhere from a few 
microseconds to a couple of seconds. Typical examples of these signals are a wrench 
falling on a metal deck or the sounds of a biological life form, such as a whale. Due to 
the diversity of these signals, there is limited a priori information that can be used to aid 
in their detection. When transients are imbedded in noise, their detection becomes very 
difficult because the energy in the noise frequently dominates the energy in the transient 
over the interval of interest. Nevertheless, in conjunction with other measurements or by 
themselves, transients can be used to estimate the track of the target submarine by 
estimating the TDOA of a transient between two sonabuoys. 

Tracking of a submarine using transients requires transient detection, 
discrimination and TDOA estimation. Prior to discrimination, a transient must be 
detected in a noisy background on the sonagram of a single sonabuoy. This thesis does 
not address the transient detection problem and therefore assumes that a transient has 
been detected on the sonogram of a sonabuoy. 

After the transient has been detected, the same transient must be found on the 
sonograms of the other sonabuoys in the sonabuoy field. This is difficult since the 
detected transient may not exist on all sonabuoys within the field. That is, the sonograms 
of the other sonabuoys could contain only noise or even a different transient in the 
window of interest. The first case (1.e., where the sonogram of the other sonabuoys 
contain only noise) is a common scenario. This happens because transients can be very 
localized and thus do not appear on the sonograms of each sonabuoy. The second case 
(1.e., where a different transient exists on the sonogram of the second sonabuoy) 1s also 
common due to the complex environment of the mission, where multiple transients from 
different sources can exist on the data records of the other sonabuoys in the field. It is 
therefore important to discriminate between different transients and to know if there is no 


transient present. If the transients arriving at two or more different sonabuoys are the 


same, then TDOA estimates can be used for target localization. In the following we first 


discuss TDOA estimation and then transient discrimination. 
B. TIME DIFFERENCE OF ARRIVAL (TDOA) ESTIMATION 


TDOA estimation is by no means simple and therefore many authors have written 
about it [Ref 1]. It 1s, however, the primary means of determining range to the target in 
passive detection, since TDOA information from multiple sonabuoys and the geometry of 
the sonabuoy field can be used to determine the position of the submarine at any 
particular instance in time. For any given value of TDOA between two sonabuoys, the 
locus of possible target positions is a hyperbolic curve. If more than one curve can be 
drawn, 1.e., if TDOA values between three or more sonabuoys can be calculated, the 
intersection of the curves determines the position of the target. 

Consider a simple two-buoy model as shown in Figure 2. The sonabuoys are 
symmetrically positioned about the origin on the x-axis as shown in Figure 2, with the 


target located at T. 





Buoy 2 (-xp,0) Buoy 1 (xp,0) 


Figure 2. Two Symmetrically Positioned Sonabuoys. 


For this geometry the time delay 7 between sonabuoy 1 and sonabuoy 2 1s given 
by 


f ty 
E =1, Eel = —— === 


C C 
(rc) - (4 (x, -x) + 2- (х, + х)? + y Js 
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where c 1s the speed of sound, D, and D, are the distances between the target and 


sonabuoy 1 and sonabuoy 2, respectively. After some algebraic manipulation’ of Eq. 2.1, 


the following expression 1s obtained [Ref 1]: 
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Equation 2.2 describes a hyperbola, and a set of hyperbolas can be drawn for different 


TDOA’s between sonabuoy 1 and 2 as shown in Figure 3. 





Figure 3. Hyperbolic curves for TDOA between two symmetrically positioned | 


sonabuoys. Sonabuoy positions are indicated by the dots. 
С. TRANSIENT DISCRIMINATION 


The objective of discrimination between transients 1s to be able to tell different 
transients apart. This is different from transient classification or transient 
characterization, which attempts to identify the source of the transient and seeks to group 


similar transients together. 


' To obtain a hyperbolic expression, Eq 2.1 must be squared, and all the cross terms must be left on the 
right hand side and all other terms on the left hand side. Both sides must now be squared again with like 
terms collected. 


Discrimination is difficult since it relies on detecting differences between the two 
received signals. These differences can be either in magnitude or phase or both. The 
problem is further complicated by the presence of noise and the fact that (due to 
differences in propagation path characteristics) the same transient arriving at two 
different locations may not look and sound the same. Achieving good discrimination at 
low SNR values is a challenging task. 

Figure 4 shows some typical cases in which transient discrimination has to take 
place. The cases shown in this figure are a localized transient such as an expanding 
bubble on a sonabuoy, a directional transient that is emitted from the hull of the 
submarine and is only received at some of the sonabuoys in the sonabuoy field, and an 


omnidirectional transient that is received at all sonabuoys in the field. 
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Figure 4. Discrimination Cases. 


Considering these cases, the TDOA cannot be estimated for all combinations of 
two sonabuoys in the field. Further, if a transient arriving at one sonabuoy 1s mistaken to 
be the same as the transient arriving at another sonabuoy, the resulting TDOA estimate 
will produce an erroneous target position. This in turn may lead to errors in tracking and 


ultimately possible erroneous SPL calculations. 


ПІ. MOMENT-BASED SIGNAL PROCESSING 


Most of the traditional work on transient processing and TDOA estimation is 
based on the use of second order statistics and classical (1.e., Fourier-based) methods of 
spectrum estimation [Ref 2], [Ref 3], [Ref 4], [Ref 5], [Ref 6]. Recently, new work has 
appeared using techniques involving higher order moments of signals and higher order 
spectra, such as the bispectrum [Ref 7]. 

The motivation for using the bispectrum in transient discrimination and TDOA 
estimation is twofold. First, the higher order spectra suppress Gaussian noise processes 
of unknown spectral characteristics. Secondly, these spectra, unlike the usual power 
density spectrum, preserve phase information [Ref 7]. In the last few years there has 
been a considerable amount of new research done in using the bispectrum and 
trispectrum for transient detection, time delay estimation and classification of signals 
Шек рз ке кегдімткеріш [Ret IIE [Ref 12]. 

This chapter defines the signal model used in this thesis and the analysis of these 
signals based on second and third moments. Finally, a set of algorithms for 


discrimination and TDOA estimation using third order moments is discussed. 
A. SIGNAL MODEL 


Before any further analysis can be done it is necessary to develop a signal model, 
which can be used for the analysis of transients arriving at two sonabuoys. As mentioned 
previously, there are three cases to be considered. In the first case, the transient arriving 
at the second sonabuoy is the same as the transient arriving at the first sonabuoy. In the 
second case, the two arriving transient signals are different while in the third case only 
noise 1s present at the second sonabuoy. It is also assumed that the sonabuoys used are 
omni-directional and that their separation in the sonabuoy field is large enough so that 
noise at the first sonabuoy is uncorrelated with noise at the second sonabuoy in both 
space and time. 

The following simple model can be used to represent a single transient arriving at 


two different sonabuoys: 


x (k)=s(k)+ m (k) 


x,(k)= As(k-L)+n,(k) Er 


where x, (k) is the noise-embedded signal arriving at sonabuoy 1, s(k) is the transient 
signal itself, and 77, (k) is white gaussian noise. Note that the signal at sonabuoy 2 is 


subject to a relative attenuation A4 and delay L with respect to the signal at sonabuoy 1. 
The frequency domain expression for these signals is 

Х.(о)- 5(ө)» М.(о) 

= (3. 2) 

X,(o)=S(oje’ + №, (о) 
where the uppercase letters represent Fourier transforms of the respective signals in Eq. 
Sli: 

If the transient arriving at sonabuoy 2 is different from the transient arriving at 

sonabuoy 1, the received signals are given by: 

2 (X) E s(k)+ 7A (2 

x, (k)=r(k-L)+m(k), 


where r(k) is the other transient arriving at sonabuoy 2 with time delay L relative to s(k). 


(3. 3) 


The corresponding expressions in the frequency domain are 
X,(@)= S(@)+ N, (o) 
| (3. 4) 
X,(@) = R(a)e™ + №, (о), 
If there is no transient present at the second sonabuoy, the two received signals 
are 
A (k) = s(k)+ n,(k) 
X) (k) mulo (k) 


In this case the signal received at sonabuoy 2 consists entirely of noise and the frequency 


(3. 5) 


domain representation of the two signals is thus: 
X, (@)= S(o)+ N,(@) 


x,()= N,(o) СА? 


В. MOMENTS AND CUMULANTS OF A RANDOM PROCESS 


For the purposes of this thesis, the following notation and definitions are used. 


The x“ moment of a real stationary random process is [Ref 7:p. 15] 


10 


meet, ..)- Elx(k)x(k +7,).x(k+7,_, )} OF 
where £ denotes statistical expectation and 7, represents the p lag. Note that m; (т) 15 


the ordinary autocorrelation function . 

In the analysis of signals using higher order statistics, cumulants rather than 
moments are generally used. Cumulants of order n are defined by certain linear 
combinations of products of moments of order 1 and lower [Ref 7:p. 15]. Cumulants of 
Gaussian random processes have the distinction that they are all zero for order n greater 
than 2. Thus signals imbedded in additive gaussian noise theoretically appear naked 
when subjected to analysis using higher order cumulants. The n” order cumulant is 
denoted by 

c ee Cum|x(k yx(k +7, Ы.Е ta N (3. 8) 
For a more detailed definition of cumulants the reader 1s referred to [Ref 7:p. 9]. 
The following relationships exist between moments and cumulants for stationary 


random processes [Ref 7:p. 9]: 
a. First Order. 
c, =m, = Eix(k)} (3. 9) 
b. Second Order. 
c; (1) = т: (т) (т! y (3. 10) 
C: Third Order 


3 3 2 2 2 3 
cT. T mam) - (mi) Im? (7,) - m; (7,) 4 m; (7, - f, )|+ 2m!) 
(3. 11) 
Using these relationships, zero mean, white Gaussian noise has the characteristics listed 


in Table 1. 





Table 1. Moments and cumulants for white gaussian noise for 7, =7, =7, =0. 


С: N-TH ORDER “MOMENTS” OF A DETERMINISTIC SIGNAL 


In this thesis, we consider transients to be deterministic signals. Consequently, 
concepts such as statistical moments are not defined. However certain operations in the 
time domain, which are analogous to estimating moments for realizations of a random 
process, are still useful for deterministic signals. Some authors (e.g., Nikias and 
Petropulu [Ref 7]) have referred to these operations as computing “moments and higher 
order spectra for deterministic signals." Since some of the techniques we have adapted 


are due to authors using this concept, we will adopt this concept here as well. 


In general the n" order moment for an energy signal, x(k), is defined as [Ref 7:p. 


78] 
т" (т.т) = У x(k)x(k +7, )..x(k+z,_,), (3. 12) 
Кш-о 
and for a power signal [Ref 7:p. 100] as 
[Ема 
mr (vsum "m >, х(К)х(&Ё+т,)....х(К+т„_,), (37195 
k=J 


where J is an arbitrary starting point of the summation and N is the period of the signal. 
From these expressions, moments can be considered as a measure of the degree of 
similarity between a signal and delayed or advanced replicas of itself. The n" order cross 


moment for n energy signals is defined as 


EA AS » XI (k)x, [XN ES ү) (3. 14) 
- ae 


The n^ order moment spectrum 1s the Fourier transform of Eq. 3.12 and can be 
expressed as Fourier transforms of the signal. For energy signals, this is given by [Ref 
7:p.85] 

M 7 (9,...9,,) 2 X(0)..X(0,,)X (o, *.. 9,4), (515) 
which can be written in terms of magnitude and phase as 

M2 (@,.-,.4)] =| X(@,)|-.| X(@,) || X(@ ++ 2,4) | 

Ф" (о »---0,-,) = G(@, )+...+ G(@,_, )— 9(0,..,,_, ), 


where M TED. ) is the magnitude term and P7 (c,...0, ,) is the phase term. 


(3. 16) 


Orders n=2, 3 and 4 are important special cases of moments. In the frequency domain 
these orders have been termed Power Spectral Density (PSD), Bispectrum, and 


Trispectrum, respectively. 
D. SECOND-ORDER MOMENTS 


1. Definitions of Moments and Spectra 


Using Eq 3.12, the autocorrelation of a deterministic signal, x(x), is written as 


m:(r)= £ х) + т) (3.17) 


=: 


From Eq 3.16, the magnitude and phase components of the PSD are given by 
M (o) - |Х(ф)` 
| (0) Kal (3. 18) 
v? (o)- 0. 


From these expressions, it can be seen that the PSD 15 an even function with no phase 


information. Similarly, the cross-correlation of two deterministic signals 1s 
т? (т)= X xk), (kc) (3. 19) 
2 Lea 
The corresponding cross-spectrum is the Fourier transform of Eq 3.19: 


Mi, (0)=X,(0)X, (0), (3. 20) 


Xx? 


which can be written in terms of magnitude and phase 


Min (о) = |х). (о) 
ж... (0) = 9, (9)- 6. (o). 


(3.21) 


2. Second Order Moments of the Received Signals 


Since cross-correlation is used extensively throughout the thesis, it is important to 
determine the cross-correlations for the three signal cases presented earlier. Before doing 


this, let us first investigate the cross-correlation between the two noise sources у, (А) апа 
17, (%) at the two sonabuoys of interest. The cross-correlation and the corresponding 


cross-spectrum are given by: 


ni, = Sn (Wnt) im 
М: (о)- зт? (с)}= N,(o)N, (o). 
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The expectation of this term 1s zero since the two white noise sources are uncorrelated. 


2 


m itself is not zero. Applying this to the case where the same 


However, the term, M 
transient arrives at the two sonabuoys, the cross spectrum 1s 


м2, (о)=|5(0) e^" -M?, (o)- Mi (v) - (3.23) 


хх; 


where M $m аге the cross terms. From this equation, it can be seen that the time delay L 


-Jet This linear phase term is therefore important 


is imbedded in the linear phase term e 
in finding the TDOA between the two signals. 
For the case where different transients arrive at the two sonabuoys, the cross- 


spectrum 1s 


M. (v) |5(0) (Је ^"^ 67 V M? (5)- M? (m) (3. 24) 


хүхэ 22 


For this case it can be seen that there is the same linear phase term e "^" as that of Eq 
3.24. Therefore the same TDOA result would be observed for Eq 3.24 and 3.25 even 
though the transients are different. The linear phase can therefore not be used as a means 
to discriminate between two transients. In the last case, where there 1s only noise present, 


the cross-correlation 1s 


M? (e)- M;, (e) (3. 25) 


771772 


For this case, it is expected that the phase will be random, and therefore no linear phase 


term will exist. 
E. THIRD-ORDER MOMENTS - BICORRELATION AND BISPECTRUM 


1. Definition of Moments and Spectra 


The third-order moment of a signal is called the bicorrelation. From Eq 3.12, the 


auto-bicorrelation of a deterministic real signal is defined as 


ту) = > x(k )x(k + 7, )x(k +7, ) (3. 26) 


k=—w 


and from Eq 3.14 the cross bicorrelation 1s 
O 2 x, (k)x, (k +7, )x,(k + 7, ). TET QN 


For this thesis, where there are only two data streams, the cross-bicorrelation will consist 


of only two signals, x,(k) and x,(k). The cross-bicorrelation functions of interest are 


m; x (г, 72 ) = E X (93 (k +T I (k +7, ) (3. 28) 
and the terms m. . x. о ub , Which are defined in an analogous way. 


The auto-bispectrum 1s 
M¿(o,,0,)=X(o, X(o,)X (e * v; ). (592) 
or in terms of magnitude and phase 


моо) -хођхе хо ко) 0000, 


0.0, )= $,(o,)+ ¢,(@,)- 9, (a, + о,) 
As can be seen, by comparing Eq 3.18 and Eq 3.30, a distinct difference between the 
second order moment spectrum and the bispectrum is that the bispectrum contains phase 


information. The cross-bispectrum is 
M?,.. (0,,0,)=X,(0,)X,(0,)X, (o, +0,) (3. 31) 


х\Х)Х3 


or in terms of magnitude and phase 


13 


Mean fessos) - боз) 





X (o, + @,) (3.32) 


T (о, ? 03 ) = ф„ (о, )+ д, (о, )- д, (о, T ©, ) 
2: Third Order Moments of the Received Signals 


Both second order moments and third order moments are used extensively 
throughout this thesis. It is therefore important to develop the third order moments for 
the three signal cases. To provide a motivation for the algorithms to be used, let us 
consider a situation in which the noise is identically zero. In reality the noise terms 
№, (о) апа №, (о) аге not zero, although the expected value of their higher order 
moments defined earlier vanish when the noise is Gaussian. 

For the case where the same transient arrives at both sonabuoys, the signals in the 
frequency domain are given by Eq. 3.2. Using Eq. 3.32, the bispectrum of these signals 


is given by 


M2... (@,,0,) =|S(@, Slo) 


ххх 


Slo, + @, ) 





(35357) 
Wiss 0:)7 s (m )- eL *ós(o.) - 6s (o. + 02), 

where ¢, is the phase of the signal S(@). Equation 3.33 is the expression for the cross- 

bispectrum between the data streams at sonabuoy 1 and sonabuoy 2 in the absence of 

noise. The signal phase terms, ¢,(@,)+¢,(@,)—¢,(@, +@,), can be eliminated by 


making use of the auto-bispectrum at sonabuoy 1, which 1s given by 


ім; (v.c. ) = IS (co, Jis(o, IS(o, +0,) 





(3. 34) 
ee (o, Tr» = 0; (о, )+ Ds (о, )- 0; (о, +@, ) 
This expression is then used to normalize the cross-bispectrum and results in the identity: 


M... (0,0, 4. 
Honor J= a ote р (5555) 
X, 172 


From Eq 3.35 it can be seen that by normalizing the bispectrum all the phase terms are 


cancelled except for the linear phase term, e "^". 


For the case of different transient arrivals, the cross-bispectrum can be obtained 


from Eq. 3.4 and Eq. 3.32 as follows 


М... (o, ? ©, ) E Ro, (о, До, T 0; ) 
нее (о, ‚©, ) = Ér (o, ) OL + ps (о, ) - Ф; (о, + 0; ) 


where ¢,(@) is the phase of the signal R(@). In this case the normalized cross- 


(3. 36) 


bispectrum takes the form 


M? , R | | 
I (o, ,0, ) = duod (o, 0,) z | (o, ) e (а) )-% (о, ) еш 


= 3. 37 
М (оа) Klo) — 


Note that, unlike the first case (see Eq. 3.35), the linear phase term cannot be separated 
from the other phase terms. In both cases however, the two-dimensional bispectrum 1s 


reduced to a function of a single variable œ. 


F. SIGNAL PROCESSING ALGORITHMS 


E Bispectrum Linear Phase Detector 


The analysis of the previous section shows that when the same signal 1s arriving 
at both sonabuoys, a linear phase term is present in the normalized function / (c, c, ). 
The time delay L can now be extracted using the following ad hoc method developed for 
time delay estimation [Ref 16]. 

To estimate the delay L, a third-order “hologram” transformation is required. 


This is defined to be [Ref 7:p. 324] 


Т(г)- T frio, o. je do do,. (3. 38) 


A TS 
Since J reduces the bispectrum to one dimension, the third-order hologram is a one- 
dimensional Fourier transform over the dimension containing the linear phase term, 
followed by an integration over the second dimension. Note that when the same signal 1s 


present at both sonabuoys, Eq. 3.38 takes the form 


Т(г)- = J Jem Pdo do. (3. 39) 


since all the elements of this second dimension are in phase they add up constructively, 
giving a strong peak at 7 - L. Since we are working in the discrete-time domain the 


third-order hologram can be rewritten as 


M -1 M -1 


T(r) ЯС И (3. 40) 


€ =00, =0 


The absolute value of 7(7) will display a strong peak at the location of the of the time 
delay between the two signals of sonabuoy 1 and sonabuoy 2. 
For the case where the two transient arrivals are different, the third-order 


hologram takes the form: 





e^ € r) ela даз do,. (3. 41) 


o= gay: 


In this case, the third order hologram contains extra phase terms dr and ds, which will 





either add in phase or out of phase. Thus, in general, T(7) will not exhibit a strong peak. 


Da Bispectrum Discriminator 


The bispectrum linear phase detector can also be used as a discriminator by 
applying a simple threshold technique to the third order hologram. The threshold 
technique is applied by first noticing, from Eq 3.41, that if the transients are different the 





hologram contains phase terms ds - ¢s - coL and magnitudes а e 1 On the other hand, if 
Slo 


the transient signals are the same, the hologram contains only the linear phase term oL 
(Eq 3.39). The phase and magnitude terms of Eq 3.41 will add constructively or 
destructively to produce peaks and valleys to the bispectral hologram. If the signals are 
the same, there will be only one peak in the hologram, and that peak will be at the delay 
Es 

In summary, one can expect that 1f the SNR 1s sufficiently large and the signals 
are the same, a peak at delay L will dominate the hologram. On the other hand 1f the 


signals are different, there will be a maximum value at delay L as well as other extrema at 
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delays other than L. These other extrema will be larger than the noise and can therefore 


be used to discriminate. A discrimination algorithm is therefore proposed as follows: 


Step | 


Estimate T(7) 


Step 2: 
Find the maximum value of the hologram 
O=ars max 7 ©) (3. 42) 
2 
(Q 1s the magnitude of the hologram at the estimated TDOA). 
Step 3: 
If z, is the value of 7 that produces the maximum in step 2 then compute 
QO; = are max{ 7 (г)- o) (3. 43) 
E 
Step 4: 
Define a threshold using Eqs. 3.42 and 3.43 as 
1 (Nae 
Q, = wale 10,00% ) (3. 44) 
Step 5: 
Compare Q, to the difference О —0О,. if 
Q; »0-O, (3. 45) 
the signals are considered to be different; if 
0,<0-0, (3. 46) 


the two signals are considered to be the same. 
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IV. SUBSPACE-BASED SIGNAL PROCESSING 


Signal subspace methods have been used in a variety of problems for estimating 
spectral lines and direction of arrival in sensor arrays. А good introduction to these 
methods for spectrum estimation can be found in [Ref 17]. Recently, the use of subspace 
methods has been explored for estimating TDOA [Ref 20], [Ref 21]. This thesis uses the 
same approach as [Ref 20] for TDOA estimation and further adapts the subspace method 
for transient discrimination. 

The discussion begins by defining the subspace model for the TDOA problem and 
then moves on to the MUSIC method, which was used for delay estimation and 
discrimination. The chapter finishes with a discussion of how the MUSIC method was 


adapted and used for discrimination 
A. SUBSPACE 


To understand the principle of subspace methods and their application to the topic 
of this thesis, consider the case where the same signal is arriving at two separate 
sonabuoys (see Eq. 3.1). The frequency domain representation of Eq. 3.1 is given by Eq 
3.2. From Eq. 3.2, it can be seen that the linear phase term e" contains the delay L 
between the transient arrivals at sonabuoy 1 and sonabuoy 2. The cross-spectrum for the 


two received signals, given by Eq 3.23, 1s repeated here for convenience 


м? (о) = |5(0) e +M7, (o)+M? (o). (4.1) 


хх 192 
It can be argued [see Ref 20] that when Eq 4.1 ıs sampled at equally-spaced values ın 
frequency, the resulting data sequence, y(k)- M 2 (w,), satisfies the conditions 
necessary to apply a signal subspace model [Ref 17]. Specifically, an NxN correlation 
matrix is estimated for the data y(k), and the corresponding N-dimensional vector space is 
divided into signal and noise subspaces. The approach followed to estimate TDOA using 


subspace methods is to project a steering vector of the form 


2] 


4(/) = |. (4. 2) 


| елемеу) | 

onto the noise subspace and plot the result as the linear phase, /, is varied. When the 
parameter / 1s equal to the true delay L between the sonabuoys, the projection onto the 
noise space is zero. In subspace algorithms, such as MUSIC, this null projection can be 
used to estimate TDOA. 

Unfortunately, the cross-spectrum for the case where different transients are 
arriving at sonabuoys | and 2 also contains a linear phase term (see Eq. 3.24). Therefore, 
the basic subspace techniques are not useful for transient discrimination without suitable 


modification. 
B. MUSIC 


In the application of subspace techniques to linear phase detection between two 
transients, the MUSIC (Multiple Signal Classification) method developed by Schmidt 
[Ref 18] is used. A brief outline of the method 1s provided here. 

The MUSIC algorithm is implemented in the frequency domain by first obtaining 


the cross-spectrum of the two signals received at the sonabuoys. Since M ae (о) 15 
formed using a DFT, this data is available at samples ®, , k= 0,1,2...Nprr -!, where Nprr 
is the size of the DFT. From this frequency domain data, a data matrix” M m is formed, 
and the correlation matrix is estimated as 


R, =M?, (m2, V. (4. 3) 


X1X%2 X1X2 
This can be decomposed into an orthonormal eigenvector matrix 


E -[E,,E (4. 4) 


me. 


and a diagonal eigenvalue matrix 


* See [Ref 17] for various methods of forming a data matrix 


22 


А,.,0 
0 ANI E >) 


The columns of Eig and E, ise form an orthogonal basis set and define the signal and 
noise subspaces, respectively. It 1s therefore possible to write R, as 


REEE E B CE A VEL (4. 6) 


sig 51р noise 
The columns of the eigenvector matrix can be used to form projection matrices for the 


signal and noise subspaces of the form [Ref 17:p. 623]: 


P, =E E? 
m O (4.7) 


When the projection matrix P is multiplied by a vector d, the result d=Pdis the 


projection of d into the corresponding subspace (see Figure 5). 


Signal 
Basis Vector 







Noise 
Basis Vector 1 





Noise 
Basis Vector 2 


Figure 5. Projection of vector d onto noise subspace by MUSIC. 


As discussed previously, subspace methods make use of the fact that the 
projection of the steering vector d(/) , given by Eq. 4.2 onto the noise subspace is zero. 


The MUSIC algorithm [Ref 17:p. 627], in particular, evaluates the quantity 


Р 1 1 
Pu == — "a OS eerat) (4. 8) 
d "(DPs dC) (Ode ed) 
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^ 


where P,, is the MUSIC indicator function for d(/)and e,...e, are the eigenvetors 


spanning the noise subspace (the eigenvector e, spans the signal-subspace which is one- 


dimensional) The value of / where P,,,, exhibits a sharp peak determines the TDOA 


between the signals. 
C. SUBSPACE DISCRIMINATOR 


A subspace discriminator was developed by comparing projections onto multiple 
subspaces using the MUSIC algorithm. The following considerations are proposed to 
motivate this discriminator. 

Since the subspace method in this thesis 1s based on the cross-spectrum between 
two sonabuoys (see Eq. 4.1) it 1s important to investigate the differences between the 
cross-spectrum for the case where the same transient arrives at two separate buoys and 
that for the case where different transients arrive at two separate buoys. This will be 


helpful in understanding the subspace discriminator. From Eq. 3.24, it can be seen that if 


the transients are different, the cross-spectrum contains phase terms ¢,()— ¢,(@)- aL 
and magnitudes ¡Rlo, Js CRE However, if the transient signals are the same, the cross- 


spectrum contains only the linear phase term wl (Eq. 3.23). Now, if the subspaces 
formed using Eq. 3.23 are compared to those formed using Eq. 3.24 the two subspaces 
would typically be rotated with respect to each other as illustrated in Figure 6. The 
implication of the rotation of the subspaces due to the phase terms of Eq. 3.24 1s that 
when the parameter / of the steering vector, d(/), is equal to the delay L between the 
sonabuoys, the projection onto the noise subspace will not be zero as in the case when the 
same transient arrives at the two sonabuoys. By comparing projections, it would be 
found that the differences in magnitude of the projections would be large if the transients 
were the same. On the other hand if the transients were different, the differences in 
magnitude between the projections would be small. These differences in magnitude are 
used to discrimnate between transients. 

This method is developed further by noticing that the noise projection matrix, 


Pose, can be written as: 
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N 
_ т "T 
SC E ER =1 ee, > (4. 9) 
{=2 


where I is an Nx N identity matrix.Let us further define P, as the projection operator for 


the subspace formed by leaving out the k” basis vector. Thus 


N 
Р, = У ее -I-e,ej. (4. 10) 


# I 
mm m 


Further we define the following set of additional indicator functions 
1 _ 1 k 
d" ()P,aQ) a^ (r-e,e; al) 


The proposed subspace-based discrimination algorithm 1s as follows: 


P, (1) = SOS (4. 11) 


5їер 1 
Estimate the correlation matrix R, and find the eigenvectors e;. 
Step 2 
Compute the maximum estimate 
Q=argmax{ ?} (4. 12) 
k,l 
Step 3 
If ko 1s the value of k that produces the maximum in step 2 then compute 
O, = аге max P) (4. 13) 
kl 
ks Kg 


Step 4 
Compare the difference O — O,. If 


О-О, <1 (4. 14) 
the signals are considered to be different. If 
Q EB О, <1 (4. 15) 


the signals are considered to be the same. 
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М. Noise 
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Noise Basis Vector 2' 
Basis Vector 2 
Figure 6. Relative rotation of subspace due to phase terms of Eq. 3.24. 
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V. SIMULATION RESULTS 


A. SIMULATION CONDITIONS 


1. Svnthetic Transients 


In this thesis a set of four synthetic transients was used to evaluate each of the 
algorithms described in the previous chapters. The four transients were a CW sinusoidal 
pulse, an exponentially decaying sinusoidal pulse, a linear phase modulated (LFM) pulse 
and a simulated finback whale transient. 

Each transient can be described as an amplitude- and phase-modulated sinusoid, 
given by the general expression [Ref 19:p. 202] 

p(t) = a(t) cos A(t), (5.1) 
where a(t) is the time-varying amplitude or “envelope” of the signal and 0(t) is the time- 
varying angle. The angle 0/1) is of the form O(t) = 2af.1 + y(t) so that 

p(t) = a(t) cos(o,t + (0), (5. 2) 
where fe is the center frequency and 11) is the phase modulation. Each of these transients 


15 described in more detail below. 


For the CW pulse, the envelope and phase modulation are given by: 


ТОЕ 
a(t) — (5959) 
0, otherwise 


y(t) 30 (5. 4) 
where T is the pulse length. A plot of this transient is given in Figure 7 for 7 = 40ms 


and f, = 50 Hz. 
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Figure 7. CW Pulse. 


The exponentially decaying sinusoidal transient is described by the following 


envelope and phase modulation terms 


here 
Do E (5. 5) 
0, otherwise 
y(t) = 0, (5. 6) 


where T is the pulse duration and ais the attenuation constant. This transient is shown in 


Figure 8 for a = 200 and T 2 40ms. 
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Figure 8. Exponentially Decaying Sinusoidal Transient. 


For the Linear Phase Modulated (LFM) transient, a Gaussian envelope of the 


form 
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y 
б О ш T (5. 7) 


0, otherwise 
was used. The phase 15 а quadratic function of time: 


Eum (5. 8) 
2T 


where Af is the desired change in instantaneous frequency over the interval 7. A plot of 


y(t) 


this transient for T = 40 ms ,o* = = and Af =100 Hz is shown in Figure 9. 
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Figure 9. LFM Pulse. 


The finback whale transient 1s the most complicated of all the transients that were 


synthesized. This transient 1s modeled as follows [Ref 8]: 


О Ота 
TB 228 


1 
alt)=! 5.9 
V) o o te 0.9) 
0, otherwise 
y(t) 2 2z - 23, e ^08: 2:Y (5. 10) 


with f. = 0. This transient is plotted in Figure 10 for T=1s. The instantaneous 


frequency decreases nonlinearly from 23 Hz to 18 Hz with the most rapid decrease 


occurring initially and the rate of decrease becoming smaller with time. 
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Figure10. | Finback Whale Transient. 
2. Signal-to-Noise Ratio 


Dunng the evaluation of the algonthms, the synthetic signals described above 
were subject to additive white Gaussian noise. The signal-to-noise ratio (SNR) is defined 
represent the average power in the signal: 


Na 
P ШТ 
transient — > x(n 


N, n=0 


as follows. Let P 


transient 


2 


(5. 11) 





Note that the transient signal power 1s normalized by the length of the transient N, and 
not the entire observation time. This 1s done to prevent the SNR from changing with 
changing observation time. The SNR in dB is then defined as 
Jm fuse 
BO NOE TOP ae Se) 
Noise noise 


where ø? is the noise variance used in the generation of the Gaussian white noise. 


noise 


B. RESULTS 


To evaluate the algonthms, two sets of experiments were used, one to perform 
TDOA estimation and the other to perform discrimination. For the first set of 
experiments, the same transient arrives at both sonabuoys. For this case all the transients 


were evaluated in turn, over a range of SNR values. For the second set of experiments, 
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all combinations of transient arrivals at sonabuoys 1 and sonabuoys 2 were tested for 
discrimination. 

In all of the experiments, the transients were set to the following lengths: 20 ms 
for the CW sinusoidal pulse, 40 ms for the exponential decaying sinusoid and LFM pulse, 
and 1 s for the whale transient. The delay, L, between the transients was always kept at 
40 samples, which 1s equivalent to 200ms. А sampling rate of four times the required 
Nyquist frequency was used for each transient. This satisfies the requirement for the 
bispectrum, which requires a sampling rate of three times the maximum frequency [Ref 
8]. The total length of each signal used was 256 samples or 1.28 s. The parameters are 


summarized in Table 2. 





PARA 
бе [5 [meme Г 


— екеш 


м [5 [memes Теё 
5 —] = | 


Table 2. Transient parameters used in experiments. 






Decaying Sinusoid 








ls TDOA Estimation 


The TDOA results will be discussed as follows. First, an example of the basic 
results using the exponentially decaying sinusoid transient will be given for both the 
subspace method and the bispectrum linear phase detector. After this, certain 
“probability” curves for TDOA will be defined and presented for each transient starting 
with the subspace method and ending with the bispectrum results. 


The subspace method uses the MUSIC algorithm with a correlation matrix (Eq. 


4.3) of size N = 60. Figure 11 shows a plot of the function Pu (see Eq. 4.8) as a 
function of / using a SNR of 12 dB. Figure 12 is a plot of the results for a SNR of 5 dB. 
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Figure 11. Subspace TDOA estimation for an exponentially decaying sinusoid 


using an SNR of 12 dB and correlation matrix size of N=60. 


As can be seen from these figures, the algorithm estimates the correct TDOA of 40 
samples or 200 ms for a SNR of 12 dB, and yields a completely incorrect value of —63 
samples for an SNR of 5 dB. 
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Figure 12. Subspace TDOA estimation for an exponentially decaying sinusoid 


using an SNR of 5 dB and correlation matrix size of N=60. 


The rationale for using a large correlation matrix is twofold. First, if the 
correlation matrix is too small, the TDOA estimate tends to be less accurate. A typical 


result is shown in Figure 13 for the 12 dB-SNR case using a correlation matrix of size N 
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— 5. Ascan be seen in this figure, a TDOA of 38 samples (instead of the correct value 40 
samples) was estimated, and the peak 1s broader than that shown in Figure 11. Secondly, 
for the subspace method to function as a discriminator, 1t was found that better results 
were achieved at lower SNR, when the correlation matrix and thus the observation space 


was larger. 
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0.8; 
057: 


| 
0.6 Estimated 
| = 
„20.5! | DTOA = 38 


0.4: | 
0.3} | 
0.2) 
0.1: 


-100 -50 0 50 100 


Figure 13. —. Subspace TDOA estimation for an exponentially decaying sinusoid 


using an SNR of 12 dB and correlation matrix size of N = 5. 


For the bispectrum linear phase detector, no windowing was applied to the data. 
Windowing is usually applied to obtain smooth estimates of the bispectrum [Ref 7:p. 
126]. Due to the short data lengths used in this thesis, however, it was found that 
windowing did not improve the results and was therefore not used. Typical results for the 
bispectrum linear phase detector are shown in Figure 14 and Figure 15 for the 
exponentially decaying sinusoid using SNR values of 12 dB and 5 dB. In both cases the 
TDOA is indicated by the maximum value of the hologram, which in tum is a single 
sample. The performance is similar to that of the subspace method. At 12 dB the result 


is exactly correct while at 5 dB the method gives completely erroneous results. 
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Figure 14. Bispectrum TDOA estimation for an exponentially decaying transient 


using an SNR of 12 dB. 
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Figure 15. Віѕресігит ТроОА estimation for an exponentially decaying transient 


using an SNR of 5 dB. 


Figure 16 shows typical TDOA results using the classical generalized cross- 
correlation methods (GCC) [Ref 2] at a SNR value of 12dB. ROTH [Ref 13], SCOT [Ref 
14] and PHAT [Ref 15] weighting was used in addition to straight cross-correlation (no 
special weighting was used). As can be seen in Figure 16 the correlation, SCOT and 


PHAT algorithms give strong peaks at the estimated TDOA while ROTH gives a very 


noisy signal with no dominant strong peak. In general, it was found that these techniques 
gave less accurate TDOA estimates than the bispectral and subspace algorithms, even at 


high SNRs. 
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Figure 16.  TDOA estimates using (a) cross correlation, (b) ROTH, (c) SCOT and 
(d) PHAT with an SNR of 12 dB. 


TDOA-estimate performance of the subspace and bispectral algorithms was 
evaluated using simulations of one thousand realizations at a fixed SNR. For each 
realization, a different white noise source was added to the original signal. The output of 
these simulations was used to create a statistical measure of performance for TDOA, 
which we call the probability of correct TDOA Py. Pr represents the fraction of 1000 
trials for which the algorithm estimates the correct time delay between transients received 
at the two sonabuoys. Pr curves were generated for SNR values of 8 dB to 20 dB for 
both the subspace method and the bispectrum linear phase detector. 

Figure 17 and Figure 18 (a) shows the Pr results of the subspace algorithm for the 
CW pulse. As can be seen from the figures, the subspace algorithm gives excellent 


results for SNR > 17 dB, where a Pr of 1.0 is achieved. 
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Figure17. Py versus SNR for CW pulse using the subspace linear phase detector. 


The Py curves for the exponentially decaying transient and the LFM transient are 
shown in Figure 18 (b) and (c). Once again, a Py value of 1.0 is achieved for SNR > 
17dB. The Pr for the whale transient (Figure 18 (d)) shows poorer results, achieving a 
maximum Py of 1.0 at a SNR of 20 dB. This poorer performance of the algorithm may 
be be attributed to two factors. First, the whale transient was the longest of all the 
transients used (being 1 s long) thus filling most of the data observation window. 
Secondly, it is a non-linear function giving it a complicated phase structure. It is also 
important to note that the results in Figure 18 (d) are based on a longer data observation 
length than used for the other TDOA results (see Table 3). 

Figure 18 (a) and Figure 19 show the Pr for the CW pulse using the bispectrum 
linear phase detector. From these figures, ıt can be seen that a Prof 1.0 ıs never achieved 
but that a Pr of 0.9 1s achieved at a SNR of 15 dB. The remaining plots in Figure 18 (b), 
(c) and (d) also show that a Pr of 1.0 is not achieved for the bispectrum method. For the 
exponentially decaying sinusoid (Figure 18 (b)) a Pr of 0.9 is achieved at 11 dB while for 
the LFM pulse (see Figure 18 (c)) a Py of 0.9 is achieved at 12 dB. As in the case of the 
subspace methods, the Pr of the whale transient was obtained using a longer data 
observation window, where a maximum Py of 0.9 was achieved at an SNR of 16 dB 


(Figure 18(d) ). 
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Figure 18. | Combined Py for subspace and bispectrum linear phase detectors: (a) 
CW pulse, (b) exponential decaying sinusoidal transient, (c) LFM pulse and (d) 


whale transient. 
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Figure 19. . Pr for the CW pulse using the bispectrum linear phase detector. 
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Figure 18 gives a summary of the TDOA results for all of the experiments. As 
can be seen in Figure 18, the bispectrum linear phase detector has higher P1 values at low 
SNR than the subspace linear phase detector for the exponential, LFM and whale 
transients. However at higher values of SNR, the subspace algorithm gives better Pr 
results (1.e., Pr of 1.0) which is never achieved using the bispectrum linear phase 
detector. 

In comparing the two methods (subspace and bispectrum) one may note that 
although the bispectrum method never achieves a Py of 1.0, in three out of the four cases 
it reaches the value of P1 = 0.9 earlier than the subspace method. Table 3 compares these 
results by listing the minimum SNR values for which a Pr of 0.9 is achieved in each of 


the test cases. 


CW Pulse LFM Pulse 


Subspace [Te 
[ Bispeetrum |IdB пав [в ев 





Table 3. Minimum SNR required to achieve Py = 0.9. 


D. Discrimination Experiments 


The transient discrimination results will be discussed in a fashion similar to that 
of the TDOA results. The discussion will therefore start by giving a brief introduction of 
how the algorithms were implemented by using the CW and LFM pulses as examples. 
After this, the probability curves for discrimination will be discussed for every 
combination of transients. 


As discussed previously, the subspace discriminator makes use of projections into 


multiple subspaces and compares the maximum values of the indicator functions P, Eq. 


4.11 with the maximum value of the function P,,, (Eq. 4.8). It is found that if the 


transients arriving at the sonabuoys are the same the difference in peak values will be 
large as shown in Figure 20(a). If the difference is found to be small (1.e., « 1) then the 
two transients are different, as shown in Figure 20(b). Thus in making these 


comparisons, it is possible to discriminate between transients. 
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Figure 20. Typical P, estimates, (a) when the signals are the same and (b) when 


the signals are different. 


In the case of the bispectrum discriminator, the phase terms ør - øs (see Eq 3.37), 
which are only present in the case of different transient arrivals, add up constructively or 
destructively thus giving features to the hologram of the bispectral linear phase detector. 
For the case where the transients are the same, the phase terms ds - ds are zero; therefore 
no extra features are added to the hologram. If the difference in magnitude between these 
features 1s large compared to the average of the hologram, it is possible to discriminate 
between transients. These characteristics can be clearly seen in Figure 21. Figure 21(a) 
shows the single peak that appears for similar transient arrivals while Figure 21(b) shows 
the multiple peaks and extra features that are present due to the phase terms dr and ds for 


different transıent arrıvals. 
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Figure 21. Bispectrum linear phase detector. Stem plot of IT (г) for (a) similar 


transient arrivals and (b) different transient arrivals. 


The experiments for transient discrimination were as follows. Each algorithm 
was evaluated using simulations of one thousand realizations at a fixed SNR. For each 
realization, a different white noise source was added to the original signal. The output of 
these simulations was used to estimate the probability of correct discrimination (Pp;). Pp, 
represents the fraction of 1000 trials for which the algorithm discriminates correctly 
between two transients. Pp; curves were generated for SNR values of 10 dB to 20 dB for 
both the subspace and bispectrum discriminators. The probability of incorrect 
discrimination, Pfa, was computed in a similar manner and also used as a performance 
measure. 

Figure 22 shows the results using the subspace discriminator where the 
exponentially decaying sinusoid is present at sonabuoy 1 and various transients are 
present at sonabuoy 2. For these curves, the vertical axis is labeled “Probability” because 


the Ps, and Pp; are shown on the same plot. In Figure 22, the Pp; curve for the exponential 
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Figure 22. Pp; and Py, (using subspace methods) between exponential decaying 
transient and exponentially decaying transient, CW Pulse, noise, LFM pulse and 


whale transient. 


decaying transient arriving at the first and second sonabuoy is shown. It can be seen 
from this curve that the algorithm discriminates correctly with a Pp; of 0.9 or larger for 


SNR values greater than 17 dB. Figure 22 also shows the Pg for the other transient 
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arrivals at the second sonabuoy (1.e., CW pulse, LFM pulse, or the whale transient) and 
the Ps, when only noise (no transient) is present at the second sonabuoy. These curves 
show that the Pg is very low. These are desirable results since they indicate that the 
algorithm can differentiate between the transients with a high probability of correct 
discrimination and a low probability of false alarm. Results for the other transients are 


summarized in Table 4; the corresponding curves are given in appendix A. It is 





“ 






important to note that the results for the whale transient are based on a longer data 
observation interval than the others. 

Transient 

Arriving 


j CW Pulse | Exponential | LFM Pulse Whale Noise 
ane 
| BUOYI | 


CW Pulse Pp; 20.35 for | Pf < 0.08 Ра <0.01 = Pr, <0.01 
EE dul ы 
Exponential | Pa < 0.02 BR >09 Р < 0.002 =0 Pa < 0.021 
тете e ee e em 








ОЕ 





t2 









LFM Pulse Pa < 0.03 Pg <0.08 Pp; 2 0.8 Bass 0] Р <0.04 
SNR>18dB 
Whale Pa <0.2 P5 «0.18 Pe < 0.18 Ры>0.8 for | Pf, <0.19 
I ol ome [o 





Table 4. Pp; results using subspace methods. 


For the bispectral discriminator, the threshold O7 (see Eq. 3.44) was adjusted by 
using a threshold gain to improve the results. The gain was applied to Eq. 3.44 in the 
following way: 


2200-0 23 (5.13) 


where G 1s the threshold gain. Threshold gains of 1, 2, 4 and 6 were used. 
Figure 23 shows the discrimination results for the exponentially decaying 
transient. using the bispectral discriminator. — Again, the vertical axis is labeled 


“Probability” because both the Pg and Pp; are shown on the same plot. For Figure 23, the 
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transient arriving at sonabuoy 1 1s the exponentially decaying sinusoid. The Pp; for this 


transient is shown in Figure 23. The other curves of Figure 23 show the results of the P, 


for the CW pulse, noise, LFM pulse and whale transient respectively. The symbols 


shown in Figure 23 are tied to the transient type and are consistant with what 1s used in 


Appendix B. The subplots in Figure 23 show the results when different threshold gains 


are used. 
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Pp; (using high-order methods) between exponentially decaying 


sinusoid transient, arriving at sonabuoy 1 and the CW pulse, exponentially decaying 


Sinusoid, noise, LFM pulse and whale transients arriving at sonabuoy 2. Threshold 


gain values of (a) 1, (b) 2, (c) 4, (d) 6 were used. 
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As can be seen in Figure 23, both Pg and the Pp; decrease with increasing 


threshold gain. A design compromise must therefore be made to have an acceptably high 
Pp; and an acceptably low Pf. For this thesis, a threshold gain of 4 was chosen. Table 5 


shows the discrimination results for the bispectral discriminator using this value for the 
threshold gain. 


Transient | P | CW Pulse | Exponential | LFM Pulse Whale Noise 
U 

Arriving | o 
Y 

| BUoY! — | 


CW Pulse 208 fori Pa <02 Pa <0.01 =0 Pg «0.001 
ai M 

Exponential Ры <0.3 Pp; 20.98 for | Р «0.15 Pj «0.001 
e eere een 

LFM Pulse | Pf <0.05 Pg «0.03 Po 2095 for Eu 0 Pg «0.001 
A Do A 

Whale Ра - 0 Paz Бе 9 Po; >0.2 for | Pf,=0 

me т. 


Table 5. Pp; using Bispectrum using a threshold gain of 4. 















The poor discrimination performance of the whale transient given in Table 5 may 
be attributed to the relatively short observation time of the data. Better results may have 
been obtained 1f longer observation times were used; however, it was not practical to run 


a large number of simulations for these longer observation times. 


С SUMMARY OF RESULTS 


The TDOA and discriminator results are summarized in Table 6 including the 
advantages and disadvantages of both algorithms. As has been previously shown in 
Figure 18, the bispectral linear phase detector tends to give a larger number of correct 
TDOA estimates at low SNR. On the other hand, the subspace algorithm gives 


consistently correct results (Pr = 1.0) at high SNR, which is never achieved by the 
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bispectrum algorithm for SNR іп the range of 8 dB to 20 dB. Therefore the more 
desirable result may depend on the application where the algorithms will be used. 

For discrimination, the bispectral method gave the highest Pp; results when a 
threshold gain of 4 was used. However, the Pf, increased with increasing SNR (see 
Figure 23). This 1s an undesirable result since even transients with high SNR will not be 
able to be separated. On the other hand, the subspace discriminator gave lower Pp; values 


but had the advantage that the Pg, remained constant for all SNR values tested. 


Algorithm | TDOA 


Bispectrum | Best Result: SNR > 11dB 
gives a Pr of 09 
(exponential transıent) 
Disadvantage: A Pr of 1 is 
never achieved at SNRs < 
20dB 

Advantage: High Prs are 
achieved at low SNRs 
Best Result: SNR > 17dB 
gives a Pr of 1. (exponential 
transient) 

Disadvantage: High Prs 
are only achieved at high 
SNRs 

Advantage: Prs of 1l are 
achieved. 





Discrimination 


Best Result: SNR > 14dB gives a 
Ppi of 0.98 (exponential transient). 
Disadvantages: Results depend on 
choice of threshold gain. Pr 
increases with increasing SNR. 
Advantages: Can achieve a high 
Pp; at relatively low SNRs. Has a 
low Pg for noise 
Best Result: SNR > 17dB gives a 
Pp; of 0.9 (exponential transient). 
Disadvantages: Computationally 
intensive. Low Pp; at high SNRs 
Advantages: Low Pg for all 
transients 











































Subspace 










Table 6. Summary of TDOA and discrimination results. 
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VI. CONCLUSIONS 


A. THESIS SUMMARY 


In this thesis, we have developed and compared two algorithms, namely the 
bispectrum and subspace linear phase detectors. These algorithms were developed for 
the purposes of transient discrimination and TDOA estimation, in order to be used as part 
of a transient tool suite to aid in the estimation of a submarine's position. Chapters III 
and IV of this thesis provide an analysis of these two algorithms. Two performance 
measures were used to evaluate the algorithms, namely the probability of correct TDOA 
(P+) and the probability of correct discrimination (Pp). 

Chapter V discusses the simulations and transients used in the evaluation of the 
algorithms. Of the four transients used, it was found that the whale transient, which is the 
longest transient, gave the worst performance. This may be due to the relatively short 
data observation times used in evaluating the algorithms. It was found that longer 
observation times produced better results for this transient; however, it was impractical to 
run a large number of cases at the longer observation times. 

It was found that both algorithms had advantages and disadvantages as 
summarized in Table 6. In general, it can be said that for TDOA the bispectral linear 
phase detector gave better results at low SNR while the subspace linear phase detector 
worked better at the higher SNR. 

For discrimination, it was found that the bispectral discriminator gave higher Pp; 
than the subspace discriminator. However, the Pg of the bispectral discriminator 
increased at higher SNR while the subspace discriminator gave a constant Pg for all SNR 
values tested. For discrimination, the advantage of having a constant Py, is desirable. 
Therefore the subspace discriminator is the best option even although it produced lower 
Pp; than the bispectral discriminator. As discussed in Chapter V, there are design trade 
offs between processing speed and performance that need to be made. For the bispectral 
linear phase detector, this trade off is in terms of threshold gain; for the subspace linear 


phase detector, this trade off is in terms of correlation matrix size. 
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B. FUTURE WORK 


The results of this thesis lead to some interesting topics for future research. These 
topics include: 
(1) The effects of the environment on a transient. In the evaluation of the two 
algorithms conducted in this thesis, no work was done to determine the sensitivity of the 
algorithms to changes in transient shape due to the environment. To do this, a 
propagation model must be used. 
(2) The effects of colored noise. In the simulations additive white Gaussian noise 
was used. This is a poor reflection of reality since ocean noise is frequently colored [Ref 
22]. It will therefore be important to see how these algorithms perform when additive 
colored noise 15 used. 
(3) Operational testing. The algorithms were tested using synthetic data. 
Ultimately, there is a need to test them on real data to obtain a measure of their 


performance and applicability in the field. 
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0.9 в 


0.4- 


0.3 — 





Figure 24. 


^ CW Pulse,Pp, | 
І 


= EXPO Decay,P;, 
Noise,P, 3 

м. ШЕМ pulse,P,_ 

E Whale,P.. 


SNR (dB) 


Pp; for the CW transient arriving at the first sonabuoy. 
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Pp; for the LFM transient arriving at the first sonabuoy. 


50 


Probability 


0.9 — | 









= CW pulse, P 
—- Whale, Ры | 
. Expo Decay, P | 

x Noise, Ee 


fa 


— | 
^^ 

| . Lfm pulse, P | 
0.3 [г “ р ы Га = 
Da2 Q ES 

ER Nm ^ A — See Ara 

ү TAN РА - ÁN Ne oyna — N aS M S 

AE dl ug err et mU cU e e усш = SE 

V =a : x — жыр ч — “ --- ES 7 = 

Kara 1; ^ Y SS е um . - 

1 `. e AD 
0.1 — + 


12 13 14 15 16 17 18 19 20 
SNR (dB) 


Figure 26. — Pp;for the whale transient arriving at the first sonabuoy. 
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APPENDIX B Pp AND Pf, USING BISPECTRUM METHOD 
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Figure 27. Pp; (using bispectral methods) between CW pulse transient , arriving 
at sonabuoy 1 and the CW pulse, exponentially decaying sinusoid, noise, LFM pulse 
and whale transients arriving at sonabuoy 2. Threshold gain values of (a) 1, (b) 2, 


(c) 4, (d) 6 were used. 
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Figure 28. Pp; (using bispectral methods) between LFM pulse transient , arriving 
at sonabuoy 1 and the CW pulse, exponentially decaying sinusoid, noise, LFM pulse 
and whale transients arriving at sonabuoy 2. Threshold gain values of (a) 1, (b) 2, 


(c) 4, (d) 6 were used. 
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Figure 29. — Py, (using bispectral methods) between whale transient , arriving at 
sonabuoy 1 and the CW pulse, exponentially decaying sinusoid, noise, LFM pulse 
and whale transients arriving at sonabuoy 2. Threshold gain values of (a) 1, (b) 2, 


(c) 4, (d) 6 were used. 
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